LISA Parameter Estimation Using Numerical 
Merger Waveforms 



g ■ J.I. Thorpe^, S.T. McWilliams\ B.J. Ke\ly\ R.P. FaheyS K. 

Q ■ Arnaud^'^, and J.G. Baker^ 

^ ■ 1 NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD 20771, 

' ^ CRESST and Department of Astronomy, University of Maryland, College Park, 

^ '. MD 20742, USA 

^ E-mail: James.I.Thorpe@nasa.gov 

^ : 

rv. 

I , Abstract. Recent advances in numerical relativity provide a detailed description 

O I of the waveforms of coalescing massive black hole binaries (MBHBs), expected to 

' be the strongest detectable LISA sources. We present a preliminary study of LISA's 

, sensitivity to MBHB parameters using a hybrid numerical/analytic waveform for equal- 

mass, non-spinning holes. The Synthetic LISA software package is used to simulate 
■ the instrument response and the Fisher information matrix method is used to estimate 

1^ I errors in the parameters. Initial results indicate that inclusion of the merger signal 

can significantly improve the precision of some parameter estimates. For example, 
. the median parameter errors for an ensemble of systems with total redshifted mass 

[ of 10^ M0 at a redshift of z ^ 1 were found to decrease by a factor of slightly more 

than two for signals with merger as compared to signals truncated at the Schwarzchild 
ISCO. 



oo 

o 



^ ■ 1. Introduction 



The inspiral and merger of massive black hole binaries (MBHBs) are expected to produce 
the strongest gravitational wave (GW) signals in the LISA frequency band. The large 
signal-to-noise ratios (SNRs) of these signals will allow LISA to not only detect, but 
measure physical parameters of the source systems. This capability will allow LISA 
to probe a wealth of astrophysical questions [1] including testing models of hierarchical 
structure formation through MBHB merger trees [2] and constraining cosmology through 
direct measurements of the distance-redshift relationship at high redshifts[3|. 

Estimates of the precision with which LISA can extract waveform parameters have 
been developing for the past several years, in parallel with progress in describing MBHB 
waveforms. Initial studies focused on equal-mass, non-spinning holesjl], while more 
recent analyses have added the effects of spin[5], E] and higher harmonics[71 [H]. All of 
these works focus on the inspiral phase of the MBHB waveform, truncating the signals 
at some time before the final merger and ringdown occur. 
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This truncation is justified given the mechanism by which the bulk of LISA's 
parameter sensitivity is expected to arise, namely modulations imparted on the signal 
from the orbital motion of the LISA constellation Although the merger signal 
contributes a significant fraction of SNR, it does so in such a brief period of time that the 
orbital modulations have little effect on the signal, reducing its usefulness for parameter 
estimation. The other, more practical reason why the merger has been excluded is that 
the post-Newtonian (PN) approximations used to compute the waveforms break down 
at or near the merger and are not expected to give accurate descriptions of the merger 
waveform. 

Numerical relativity (NR) has now provided a description of the complete waveform 
from inspiral through merger and ringdown[9l [10], at least for a restricted set of 
parameters. Analyses using hybrid PN/NR waveforms confirm that the merger portion 
of the signal often dominates the total SNR of a LISA MBHB observation [llj. These 
waveforms make it possible to check the assumption that the merger has a minimal 
effect on parameter est imat ion [T2l [13] . In this paper we present an initial study of LISA 
parameter estimation using complete waveforms for equal-mass, non-spinning MBHBs. 

2. Methodology 

2.1. Waveform 

The waveform used in these parameter estimation studies is a hybrid of PN and NR 
waveforms as described in section IV of Baker, et al .[IT]. Briefly, the PN waveform 
(3.5PN in phase and 2.5PN in amplitude) spans the time from 2 x 10^ Mrest{G / c^) 
to 328Mrest{G / c^) before the merger, where Mrest is the total mass of the MBHB. At 
328Mrest{G / c^) , the waveform is stitched to a NR simulation that extends through 
merger and ringdown. The matching time is chosen at the point where the estimated 
accuracy of the NR simulations exceed those of the PN expressions. 

Since the hybrid waveform relies on NR simulations to generate the merger segment, 
it is not practical to use this technique to generate the large sets of waveforms with 
varying parameters that are needed for many data analysis techniques. However, if we 
restrict ourselves to variations in extrinsic parameters and keep intrinsic parameters 
such as mass ratio and hole spins fixed, it is possible to use a single "master" waveform 
and vary the extrinsic parameters using analytic formulae. 

For our purposes, the desired waveforms are dimensionless strain in the Solar 
System Barycenter (SSB) frame, hssB- These depend on six extrinsic parameters: the 
redshifted total system mass M = Mrest{^ + z), luminosity distance D^, coalescence time 
tc, and three angles describing the orientation of the binary, in this case an inclination 
L, initial orbital phase 0o, and a polarization ip. It is convenient to construct a rescaled 
complex strain in the source frame, hgource = Rext {h+ + ihx), which scales with M^est, 
has a constant amplitude in the far-field limit, and can accommodate the necessary 
angular dependencies through rotations in the complex plane. Rext is a conventional 
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notation for the radius at which the radiation has been extracted in the numerical 
simulation, although in our case, it can equally well be used as a point of reference for 
rescaling the PN amplitude. With these definitions, we can use the following relationship 
between the rescaled strain in the source frame and the dimensionless strain in the SSB, 
hssB- 

^5 G M 




^SSB — \ oTT^ 



hsourceCOs' 6^'^° + /^L^, siu^ (0 6"^^ 



Two additional extrinsic parameters describe the sky location of the binary in the 
SSB frame. In this case we use galactic longitude $ and latitude 9. The effect of these 
parameters of the waveform is handled by the instrument model, described in the next 
section. In total, our model uses eight extrinsic parameters to determine hsss from 
hsource- Ih the remainder of this work, the vector A = (M, D^, 0, $, l, 0o, ip, t^) denotes 
these parameters and the shorthand notation h = h{X) represents the transformation in 

m . 



2.2. Instrument Model 

There are two necessary components for modeling the sensitivity of a gravitational wave 
detector [H]: a response function relating the gravitational wave signal to the instrument 
output and a noise model describing all other components of the instrument output. 

In the case of LISA, the instrument outputs are the Time Delay Interferometry[T5l 
[16] (TDI) streams. TDI is a technique by which the differential phase measurements 
made at each spacecraft are combined with appropriate time delays in order to form 
outputs for which laser phase noise is canceled but GW signals remain. For data 
analysis purposes, it is convenient to use a set of "optimal" TDI channels |17j that 
form an orthogonal basis in the TDI signal space. For an optimal set, each channel 
contains independent noise and can therefore be treated as a separate detector. We 
have defined a set of optimal channels A = {A' , E' ,T') where A' = — X), 

E' = + X - 2 ■ F), and r = \^{Z + X + Y) where (X, F, Z) are the first- 

generation Michelson TDI variables. The primed notation is used to distinguish A from 
the {A, E, T) used by Prince, et al . [17] . 

The Synthetic LISA software package [T5] is used to model the instrument response 
to the GW signal. Synthetic LISA models the orbits of the individual spacecraft 
and computes the constellation geometry. The GW response in each arm is then 
computed and assembled to form the appropriate TDI combinations. Although the 
computational cost of a numerical simulation exceeds that of an analytic model, the 
Synthetic LISA model provides higher fidelity, particularly at high Fourier frequencies 
where the instrument response becomes complex and time-dependent. 

It would also be possible to use Synthetic LISA to compute the noise output in 
the A channels as well as the signal response. The difficulty with this approach is that 
statistical variations cause the noise curve to change from one run to the next, washing 
out the subtle changes corresponding to variations in waveform parameters. This effect 
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could be mitigated by averaging repeated simulations with varying seeds, however this 
would add to the computational cost. 

An alternative approach, the one taken by most estimates of LISA sensitivity in 
the literature, is to derive an analytic expression for the average noise level in a TDI 
channel given the noise spectra of the individual components. For A, the one-sided 
power-spectral densities (PSDs) of the noise are given by 

Sa',e' = 2 sin^ (C) ■ {2 ■ [3 + 2 cos (C) + cos (2C)] 

+ [2 + cos(C)]5,,(/)}, (2) 
St' =8sin2(C)-sin2(C/2)- [4sin2(C/2)^p^(/) + 5op(/)], (3) 

where ( = 27r/L/c, / is the Fourier frequency, L is the average detector arm- length, 
and Spmif) and Sop{f) are the strain-equivalent proof- mass acceleration and optical- 
path noises of the LISA links respectively. For this work, Spm{f) = (2.5 x 10~^^ Hz~^) • 

(//I Hz)-' ■ ^l + (//0.1mHz)-' and Sop{f) = (L8 x lO'^^ Hz"!) ■ (//I Hz)^ 

In addition to the noise introduced by the instrument, MBHB signals are also 
partially obscured by a "foreground" of other GW signals. This foreground is composed 
of a superposition of tens of thousands of compact binaries in our galaxy and its 
neighborhood. We estimate the noise contribution in A' and E' of this stochastic 
foreground as Sgai{f) = [4^sin(^)]' ■ Sconfif), where Sconfif) is taken from equation 
(14) of Timpano, Rubbo, and Cornish[19j. The pre-factor in Sgai{f) is an estimate of 
the signal response functions of A' and E' for the frequencies where Sconfif) is defined, 
40 /^Hz < / < 8 mHz. At these low frequencies, the response in the T' channel is near 
zero. Consequently no foreground is added to the T' channel. 

One downside to using different instrument models for the signal and the noise is 
that artifacts can arise when these models are inconsistent. For example, the analytic 
formulas predict zero noise at /„ = nc/L, n = 1,2,3.... Were analogous formulae used 
to predict the signal response, the result would be zero signal at / = /„. In practice, 
the Synthetic LISA model predicts some finite signal at / = /„, partly due to errors 
in estimating PSDs from discrete time series. The result is an infinite contribution 
to the SNR at / = /„. To alleviate this problem, we added an artificial "floor" of 
(10-^° Hz"^) ■ (//lHz)2 to the noise models for A. 

The T' channel has a related problem at low frequencies, where St' oc The 
corresponding signal response function should also fall off with at low frequencies in 
order for SNR to converge. In practice, the T' channel response for our signals leveled 
off at a finite value, leading to an erroneous increase in SNR (and parameter sensitivity). 
It may be possible to address the problem with the T' channel by adding an artificial 
noise floor to Sr'if) or by only using it in the high-frequency limit. For this work, 
however, we have chosen to neglect the contributions of the T' channel until the issue 
can be further investigated. 

Finally, we note that these types of artifacts may arise when actual instrumental 
data is analyzed. Most analysis methods require some model of the instrument response 
and noise, each of which will have some discrepancies with the actual instrument 
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behavior. 

2.3. Parameter Sensitivity Calculation 

With the waveform described in Sec. 12.11 and the instrument model described in Sec. 
I2.2l in hand, the next step is to develop a procedure for evaluating parameter sensitivity. 
Like many others in the community, we have selected the Fisher information matrix 
approach [20]. 

The steps for the approach are as follows. A SSB waveform with a particular set 
of extrinsic parameters Ao is generated from the master waveform, ho = hssB {^oj- 
A set of perturbed waveforms is then generated by perturbing each of the parameters 
individually, hi = hsssiK + ^Ki), where SXi is the size of the perturbation in the i 
parameter direction. The instrument response to ho and hi is computed using Synthetic 
LISA. The partial derivative of the instrument response with respect to each parameter 
is then estimated using a one-sided difference, 
dA 

^ A=A. 

where A,, is the response to ho and Aj is the response to hi. These parameter derivatives 
are used to estimate the Fisher information matrix 



Ai Ao 



/dA 



OA 



(5) 



where (a|6) denotes the noise-weighted inner product, 

H.).2f«±ip^^ (6) 

The quantity Sn{f) is the PSD of the appropriate TDI channel. In our case we estimate 
the Fourier transform of dA/dX from the time series using a discrete Fourier transform 
with tapering applied to the early and late portion of the signal to reduce spectral 
artifacts. The inner product in ([6]) is evaluated as a discrete sum with finite upper and 
lower integration limits. 

Since the A channels are constructed to be linearly independent, it is appropriate 
to sum their individual Fisher matrices to yield the effective Fisher matrix for the 
instrument as a whole. However, as discussed in section 12. 2^ the T' channel was not 
well-behaved in our instrument model. For that reason, we have omitted the T' channel 
and defined the composite Fisher matrix as Tij = T^^ ^ + F^-^ \ Taking the usual limit 
of large SNR, we compute the covariance matrix E*-' ^ (F^^)*'' and parameter errors 
AA* = VT/^. In order to explore the effect of the choice of Aq, we perform multiple runs 
with varying Ao and estimate the probability density functions of S*-' and AA*. 

3. Results 

For our initial trials, we choose an ensemble of 140 systems, each with M = 10^ Mq and 
Dl = 6 Gpc, corresponding to z ~ 1 in standard A-CDM cosmology. The remaining 
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Table 1. Mean and median values for parameter errors and SNR for an ensemble 
of 140 equal-mass, non-spinning binaries with AI = 10® at Dl — 6Gpc {z ~ 1). 
For 'pre-ISCO' , integrations of Fisher matrix elements were cut off at a frequency 
corresponding to Schwarzchild ISCO. 





total mean 


pre-ISCO mean 


total median 


pre-ISCO median 


SNRae 


1.12 X 10^ 


8.35 X 10^ 


9.28 X 10^ 


6.98 X 10^ 


Aln(M/M0) 


2.35 X 10-6 


4.00 X 10-6 


1.66 X 10-6 


3.26 X 10-6 


Aln(DL/lpc) 


3.78 X 10-1 


6.07 X 10-1 


4.48 X 10-2 


1.18 X 10-1 


Ae (deg) 


3.84 


7.19 


2.15 


5.28 


A$ (deg) 


17.7 


26.9 


4.64 


9.19 


At (deg) 


81.9 


137 


1.08 


3.12 


A0O (deg) 


607 


1270 


1.59 


4.26 


Alp (deg) 


620 


1290 


7.46 


15.9 


Atc(sec) 


17.1 


39.5 


6.87 


27.3 



parameters (B, $, l, (po, "^P, and t^) were randomly varied over the ensemble. For each 
member of the ensemble, Fjj was computed using two different upper frequency limits 
for the inner product in ([6]). In one case, the integral was truncated at the Schwarzchild 
inner-most stable circular orbit (ISCO) frequency whereas in the second case it was 
continued for an order of magnitude beyond the peak (merger) frequency. Table [1] 
lists the mean and median values of the SNR and parameter errors for the ensemble. 
Histograms of each parameter error for both cases are shown in Figure [H Note that the 
plotted quantities are the direct results of the technique outhned in [21 without imposing 
any physical limits on the sizes of the errors. For example, an error in galactic latitude of 
~ 1000 deg should be interpreted as the measurement providing no useful O information 
for that system. 

The median improvement in A A' upon inclusion of the merger is a factor of ~ 2, with 
the exception of At^, which improves by a factor of ~ 4. In all cases, the improvement 
factor is greater than that for the median SNR of ~25%, suggesting that most of the 
improvement is via a mechanism other than direct increase in SNR. 

One possible mechanism is that a pair of parameters that is nearly degenerate in the 
inspiral phase becomes decoupled in the more complex merger phase of the waveform. 
Table [2] lists the percentage decrease in covariance for each parameter pair from pre- 
ISCO to total waveform. Note that all of the covariances decrease when the post-ISCO 
waveform is added. Covariances with the coalescence time exhibit the most dramatic 
decrease. This can be qualitatively explained by arguing that the merger provides a 
sharp feature that constrains the coalescence time more tightly and breaks degeneracies 
with other parameters. 
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Figure 1. Histograms of SNR and AA' for ensemble summarized in Table [TJ The 
dashed lines are the results with a cutoff at ISCO, the solid lines include the entire 
waveform. 



Table 2. Median percentage reduction in covariance from pre-ISCO to total signal, 
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83.5 


83.4 


95.4 


e 






80.9 


89.6 


80.8 


77.2 


81.9 


97.6 


$ 








82.9 


83.1 


80.8 


81.8 


87.2 


L 










76.8 


87.0 


86.6 


94.8 














78.7 


78.5 


91.8 
















78.3 


93.0 


tc 
















94.2 



LISA Parameter Estimation Using Numerical Merger Waveforms 



8 



3.1. Evolution of Parameter Errors 

The results in Table [11 Figure [H and Table [2] demonstrate the improvement in parameter 
estimation when the merger and ringdown portions of the waveform are included. A 
related question is to ask how the parameter information accumulates with time. This 
is of great interest for the problem of searching for electromagnetic (EM) counterparts 
for MBHB mergers. If the final portion of the signal provides significant improvement in 
parameter estimates, there will be pressure to reduce the latency between LISA and EM 
observatories. This could affect decisions on data downlink operations, inter-spacecraft 
communications, and other aspects of the mission. 

Figure [2] shows the effect of varying the upper cutoff frequency of the inner 
product in ([6]) on the SNR, AD^, AO, and A$. From the EM counterpart search 
perspective, these (along with perhaps tc) are the most critical parameters. Since our 
waveforms are composed of a single mode, the GW frequency and {t — tj can be related 
analytically [2T]. Included on the plots in Figure [2] are points corresponding to (t — tc) = 
1 min., 1 hr., 1 day, and 1 wk. Figure [1] indicates that A0 and A$ improve by more 
than an order of magnitude in the final day before merger. Lang & Hughes |22] looked 
at spinning, precessing, unequal mass binaries using PN waveforms and found that sky 
position error decreased by similar amount during the final day before merger. 

3.2. Effect of Mass 

As a test of the effect of total redshifted system mass on our parameter estimation 
results, we simulated an ensemble of 140 equal-mass, non-spinning MBHBs with 
M = 10^ M0 and Dl = 6 Gpc. Table [3] gives a summary of the mean and median 
parameter errors for this run. Comparing with the analogous ensemble of M = 10^ Mq 
systems summarized in Table [Tj the heavier systems generally exhibit a reduced SNR 
and increased parameter errors. This can be explained by the shift of the GW signal to 
lower frequencies and its obscuration by the steeply-rising noise of the LISA detector. 
As with the M = 10^ Mq systems, the reduction in parameter errors with the inclusion 
of the merger signal is still present in the M = 10^ Mq case. 

4. Conclusion 

The results presented above indicate that inclusion of the merger signal can have a 
significant impact on LISA's ability to measure source parameters. It is interesting 
to compare this effect with other effects that have been studied previously. In Table 
m we compare our results for M = 10^ M0 at 2; ^ 1 with prior published works. 
The parameters compared are the fractional luminosity distance ADl/D^ ~ Aln(I?x,) 
and the solid angle subtended by the error elhpse, AQtv. The latter quantity can be 
computed from the sky position angles, their errors, and their covariance. Aside from 
the differences in the waveform physics (spin, merger, etc.), the biggest difference in the 
studies is the length of the waveforms. Ours are a factor of ~ 30 shorter than those 
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Figure 2. Effect of upper frequency cut on SNR, AD^, AQ, and A$ for the ensemble 
summarized in Table [TJ Solid line denotes median values while shaded spans 20*'' 
to 80*'' percentile range. The vertical dashed line denotes the Schwarzchild ISCO 
frequency. Individual solid markers denote time before merger. 



Table 3. Same as Table [T] for an ensemble of 140 equal-mass, non-spinning binaries 
with M = 10^ Mq at Z^l ^ 6 Gpc (z - 1). 





total mean 


pre-ISCO mean 


total median 


pre-ISCO median 


SNRae 


98.5 


26.8 


85.2 


22.9 


Aln(M/M0) 


1.55 X 10-6 


2.66 X 10-6 


1.01 X 10-6 


1.89 X 10-6 


Aln(DL/lpc) 


2.17 


2.83 


8.04 X 10-1 


1.07 


Ae (deg) 


22.4 


33.5 


19.8 


27.6 


A$ (deg) 


58.4 


103 


22.0 


35.1 


Ai (deg) 


9.65 


25.2 


5.88 X 10-1 


1.40 


A0O (deg) 


1.58 X 10^ 


2.80 X 10^ 


33.3 


43.9 


Alp (deg) 


1.62 X 10^ 


2.87 X 10^ 


79.9 


111 


Atc(sec) 


14.7 


77.2 


13.1 


68.5 
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Table 4. Comparison of source location for equal-mass MBHBs at z ^ 1 in this work 
to prior published works. S,P, and M refer to waveforms with spin, spin precession, 
and merger respectively. (* estimated from the results in Table II of ^jj ^ estimated 
from histograms in FIG. 2 of 5J) 





M(106 Mo) 


duration (yr) 


S 


P 


M 


AQn (deg') 


ADl/Dl 


Cutler [4J 


4.0 


1.0 


no 


no 


no 


4 X 10-^* 


7 X 10-2* 


Vecchio[5j 


4.0 


1.0 


no 


no 


no 


6 X lO-^t 


2 X 10-2t 


Vecchio|5j| 


4.0 


1.0 


yes 


no 


no 


3 X lO-^t 


6 X 10-3t 


Lang & Hughes [22] 


1.2 


~1.5 (avg) 


yes 


yes 


no 


5.56 X 10-2 


3.57 X 10-3 


this work 


1.0 


0.03 


no 


no 


no 


5.08 X 10^ 


1.18 X 10-1 


this work 


1.0 


0.03 


no 


no 


yes 


1.13 X 10^ 


4.48 X 10-2 



used in the other pubhshed works due to the computational effort required to simulate 
long waveforms with Synthetic LISA. 

The best comparison that can be made with this work is to compare the non- 
spinning results of Cutler and of Vecchio to our pre-ISCO case. AD^/Dl is a factor 
of 2 ~ 5 larger in our results than the published results. When the merger is added, 
the three results for AD^/Dl approximately agree. This suggests that observing the 
final ~ 10 days of a M = 10^ Mq MBHB at 2; ~ 1 will constrain Dl equally as well as 
observing the entire year prior to ISCO. 

The results for AQn for the short waveforms are considerably worse. This is 
consistent with the picture that modulations due to LISA's orbital motion provide 
much of the sensitivity to these parameters. Nevertheless, the improvement in Afi^v 
upon inclusion of the merger is a factor of ~ 5, roughly the same factor as the inclusion 
of spin precession. 

Whether this same level of improvement will persist when long-duration 
observations, spin, and precession are added will depend on the exact mechanism by 
which the improvement arises. If both improvements result from breaking the same 
parameter degeneracy, then the improvement is likely to be minimal. If however, the 
mechanisms for improvement are orthogonal, it is possible that the relative improvement 
upon including the merger will be similar to the case here. 
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